
%Take draws from the parameter space. 
%These draws are held in BETAS which is size NPxNVxNDRAWS
%Then calculate the probability of each person's sequence of choices at
%each draw for that person. 
%These probabilities are held in PROBS which is size NPxNDRAWS

BETAS=randi(NGridPts,[NP,NV,NDRAWS]); %BETAS go from 1 to NGridPts
BETAS=(BETAS-1)./(NGridPts-1);   %Now BETAS go from zero to 1 inclusive
for r=1:NV
  BETAS(:,r,:)=COEF(r,1)+(COEF(r,2)-COEF(r,1)).*BETAS(:,r,:);  %Now BETAS go from lower limit to upper limit for each coefficient
end;

rrate = squeeze(BETAS(:,1,:));
QH = squeeze(BETAS(:,2,:));
Psi = zeros(NP,T+1,NDRAWS);
%%%% 
Psi(:,1,:) = 1;
for ii=1:T
  Psi(:,(ii+1),:)= QH .* ((1 ./(1+rrate)).^(ii));
end
    
%expand to correct number of rows using the ID column
Psi = Psi(XMAT(:,1),:,:);

% load values of discounted variables
%FirstRun will change to zero in bootstrapping, so rows align with
%bootstrap sample
if (FirstRun==1 & nset==3)
    %pricetime = csvread('price.time.csv', 1)./10;
    captime = csvread('CAP.time.csv', 1); %no difference for ana model
    swtime = csvread('SW.time.csv', 1);%no difference for ana model
    
    pricetime = csvread('price.time.anafull.csv', 1)./10;
    buffertime = csvread('buffer.time.anafull.csv', 1); %full attendence
    jobstime = csvread('jobs.time.anafull.csv', 1);
    infratime = csvread('infra.time.anafull.csv', 1);
    wlifetime = csvread('wlife.time.anafull.csv', 1);
    qualitytime = csvread('quality.time.anafull.csv', 1);
    
    ppricetime = csvread('price.time.anapart.csv', 1)./10;
    pbuffertime = csvread('buffer.time.anapart.csv', 1); %partial attendence
    pjobstime = csvread('jobs.time.anapart.csv', 1);
    pinfratime = csvread('infra.time.anapart.csv', 1);
    pwlifetime = csvread('wlife.time.anapart.csv', 1);
    pqualitytime = csvread('quality.time.anapart.csv', 1);
     
    igpricetime = csvread('price.time.anaig.csv', 1)./10;
    igbuffertime = csvread('buffer.time.anaig.csv', 1); %no attendence
    igjobstime = csvread('jobs.time.anaig.csv', 1);
    iginfratime = csvread('infra.time.anaig.csv', 1);
    igwlifetime = csvread('wlife.time.anaig.csv', 1);
    igqualitytime = csvread('quality.time.anaig.csv', 1);

elseif (FirstRun==1 & nset==2)
    captime = csvread('CAP.time.csv', 1); %no difference for ana model
    swtime = csvread('SW.time.csv', 1);%no difference for ana model
    
    pricetime = csvread('price.time.ana.csv', 1)./10;
    buffertime = csvread('buffer.time.ana.csv', 1); %full attendence
    jobstime = csvread('jobs.time.ana.csv', 1);
    infratime = csvread('infra.time.ana.csv', 1);
    wlifetime = csvread('wlife.time.ana.csv', 1);
    qualitytime = csvread('quality.time.ana.csv', 1);
    
  
    igpricetime = csvread('price.time.anaig.csv', 1)./10;
    igbuffertime = csvread('buffer.time.anaig.csv', 1); %no attendence
    igjobstime = csvread('jobs.time.anaig.csv', 1);
    iginfratime = csvread('infra.time.anaig.csv', 1);
    igwlifetime = csvread('wlife.time.anaig.csv', 1);
    igqualitytime = csvread('quality.time.anaig.csv', 1);
    
    
end

dbenefits = zeros(NROWS, 7, NDRAWS);
dprice = zeros(NROWS, NDRAWS);

%BETAS2 = BETAS(:,[9,12,15,18,21],:);
%BETAS1 = BETAS(:,[3,4,8,11,14,17,20],:);
%BETAS3 = BETAS(:, [10,13,16,19,22],:);
for ii=1:NDRAWS
        dprice(:,ii) = sum(squeeze(Psi(:,:,ii)) .* pricetime, 2);
        dbenefits(:,1,ii) = sum(squeeze(Psi(:,:,ii)) .* captime, 2);
        dbenefits(:,2,ii) = sum(squeeze(Psi(:,:,ii)) .* swtime, 2);
        dbenefits(:,3,ii) = sum(squeeze(Psi(:,:,ii)) .* buffertime, 2);
        dbenefits(:,4,ii) = sum(squeeze(Psi(:,:,ii)) .* infratime, 2);
        dbenefits(:,5,ii) = sum(squeeze(Psi(:,:,ii)) .* jobstime, 2);
        dbenefits(:,6,ii) = sum(squeeze(Psi(:,:,ii)) .* qualitytime, 2);
        dbenefits(:,7,ii) = sum(squeeze(Psi(:,:,ii)) .* wlifetime, 2);
end
%v=sum(bsxfun(@times,dbenefits,BETAS(XMAT(:,1),[3,4,8,11,14,17,20],:)),2);
%sum(NROWSx(NV-1)xNDRAWS,2) gives NROWSx1xNDRAWS
if(nset==3)
    v=sum(bsxfun(@times,dbenefits,BETAS(XMAT(:,1),[3,4,8,11,14,17,20],:)),2);
    v=squeeze(v);
    scale = BETAS(XMAT(:,1),5,:); %+ BETAS(XMAT(:,1),6,:) + BETAS(XMAT(:,1),7,:);
    v=bsxfun(@minus,v,squeeze(scale).*dprice); %Subtract price
    v=squeeze(v);
elseif(nset==2)
    %%%CHANGE INDICIES FOR SET3
    v=sum(bsxfun(@times,dbenefits,BETAS(XMAT(:,1),[5,6,7,9,11,13,15],:)),2);
    v=squeeze(v);
    scale = BETAS(XMAT(:,1),3,:); %+ BETAS(XMAT(:,1),6,:) + BETAS(XMAT(:,1),7,:);
    v=bsxfun(@minus,v,squeeze(scale).*dprice); %Subtract price
    v=squeeze(v);
end

if(nset==3)
    dbenefits = zeros(NROWS, 5, NDRAWS);%change dimension
    for ii=1:NDRAWS
        dprice(:,ii) = sum(squeeze(Psi(:,:,ii)) .* ppricetime, 2);
        dbenefits(:,1,ii) = sum(squeeze(Psi(:,:,ii)) .* pbuffertime, 2);
        dbenefits(:,2,ii) = sum(squeeze(Psi(:,:,ii)) .* pinfratime, 2);
        dbenefits(:,3,ii) = sum(squeeze(Psi(:,:,ii)) .* pjobstime, 2);
        dbenefits(:,4,ii) = sum(squeeze(Psi(:,:,ii)) .* pqualitytime, 2);
        dbenefits(:,5,ii) = sum(squeeze(Psi(:,:,ii)) .* pwlifetime, 2);
    end
    %%% CHANGE INDICES
    v3=sum(bsxfun(@times,dbenefits,BETAS(XMAT(:,1),[9,12,15,18,21],:)),2); 
    v3=squeeze(v3);
    scale = BETAS(XMAT(:,1),6,:); %+ BETAS(XMAT(:,1),6,:) + BETAS(XMAT(:,1),7,:);
    v3=bsxfun(@minus,v3,squeeze(scale).*dprice); %Subtract price
    v3=squeeze(v3);%NROWSxNDRAWS
end

dbenefits = zeros(NROWS, 5, NDRAWS);%change dimension/CLEAR
for ii=1:NDRAWS
        dprice(:,ii) = sum(squeeze(Psi(:,:,ii)) .* igpricetime, 2);
        dbenefits(:,1,ii) = sum(squeeze(Psi(:,:,ii)) .* igbuffertime, 2);
        dbenefits(:,2,ii) = sum(squeeze(Psi(:,:,ii)) .* iginfratime, 2);
        dbenefits(:,3,ii) = sum(squeeze(Psi(:,:,ii)) .* igjobstime, 2);
        dbenefits(:,4,ii) = sum(squeeze(Psi(:,:,ii)) .* igqualitytime, 2);
        dbenefits(:,5,ii) = sum(squeeze(Psi(:,:,ii)) .* igwlifetime, 2);
        
end 
    
if(nset==2)
    v2=sum(bsxfun(@times,dbenefits,BETAS(XMAT(:,1),[8,10,12,14,16],:)),2); 
    v2=squeeze(v2);   
    scale = BETAS(XMAT(:,1),4,:);%NROWSxNDRAWS
    v2=bsxfun(@minus,v2,squeeze(scale).*dprice); 
    v2=squeeze(v2); 
elseif(nset==3)
    v2=sum(bsxfun(@times,dbenefits,BETAS(XMAT(:,1),[10,13,16,19,22],:)),2); 
    v2=squeeze(v2);   
    scale = BETAS(XMAT(:,1),7,:);%NROWSxNDRAWS
    v2=bsxfun(@minus,v2,squeeze(scale).*dprice); 
    v2=squeeze(v2); 
end

%Add up the non-price variables times their WTP
v = v + v2;
if(nset==3) 
    v = v + v3;
end

v=exp(v);
v(isinf(v))= 1.8e307;
sparsematrix=bsxfun(@eq,sparse(1:NCS)',XMAT(:,2)'); %NCSxNROWS
denom=double(sparsematrix)*v;                           %NCSxNDRAWS
denom(denom==0)=0.0000000000000001;
p=v(XMAT(:,3)==1,:)./denom;                     %NCSxNDRAWS
p(p==0)=0.0000000000000001;
personid=XMAT(XMAT(:,3)==1,1);                  %NCSx1
sparsematrix=bsxfun(@eq,sparse(1:NP)',personid');    %NPxNCS
PROBS=double(sparsematrix)*log(p);                      %NPxNDRAWS
PROBS=exp(PROBS);                               %NPxNDRAWS
clear sparsematrix v v2 v3denom p personid dprice dbenefits Psi rrate scale
if WantBoot==0
    clear igpricetime igqualitytime igwlifetime iginfratime igjobstime igbuffertime igswtime igcaptime
    clear qualitytime wlifetime infratime jobstime buffertime swtime captime
    if nset==3
        clear pqualitytime pwlifetime pinfratime pjobstime pbuffertime
    end
end


